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A  folding  wing  structure  consisting  of  three  components:  fuselage,  in¬ 
board  wing  and  outboard  wing,  is  modeled  computationally  using  a  geo¬ 
metrically  nonlinear  structural  dynamics  theory  based  upon  von  Karman 
strains  and  a  three-dimensional  nonlinear  potential  flow  aerodynamic  model. 

The  structural  dynamic  equations  of  motion  are  discretized  in  space  using 
a  discrete  Ritz  basis  derived  from  finite  element  analysis  and  component 
synthesis  and  the  aerodynamic  model  is  discretized  using  a  vortex  lattice. 
Results  from  the  computational  model  are  compared  to  those  from  exper¬ 
iments  designed  and  tested  in  the  Duke  University  wind  tunnel  for  three 
folding  wing  configurations.  It  is  found  that  overall,  limit  cycle  oscilla¬ 
tion  magnitude  and  frequency  results  from  theory  compare  well  with  those 
measured  in  the  experiment.  The  study  also  indicated  that  as  the  flow  ve¬ 
locity  is  increased,  the  computed  limit  cycle  response  of  the  inboard  wing 
showed  a  higher  level  of  dynamic  complexity,  as  measured  by  the  number 
of  significant  frequencies  contained  in  the  response,  than  the  experimen¬ 
tal  model.  The  theoretical  model  also  predicts  that,  for  the  two  folding 
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wing  configurations  with  the  smaller  outboard  folding  angles,  the  outboard 
wing  limit  cycle  oscillation  tip  displacement  results  are  very  similar  and 
significantly  larger  in  magnitude  than  those  of  the  configuration  with  the 
largest  outboard  folding  angle.  Overall  both  the  theoretical  and  experi¬ 
mental  models  demonstrate  that  different  flutter  and  limit  cycle  behavior 
may  occur  for  different  folding  wing  configurations.  Also  it  appears  that 
structural  nonlinearities  are  stronger  than  aerodynamic  nonlinearities  for 
the  cases  studied. 


I.  Introduction 

Efforts  to  develop  morphing  air  vehicles  with  multiple  mission  capabilities  have  recently 
been  undertaken  by  several  research  teams  including  NASA’s  Aircraft  Morphing  program1 
and  the  Defense  Advanced  Research  Projects  Agency’s  Morphing  Aircraft  Structures  pro¬ 
gram.2  One  such  morphing  wing  structure  is  the  folding  wing  concept.  With  multiple 
individually  articulated  sections,  various  wing  geometries  can  be  achieved  in-flight  allowing 
for  multirole  missions  with  the  same  aircraft. 

Several  parametric  studies  have  been  performed  to  assess  the  linear  aeroelastic  charac¬ 
teristics  of  generic  folding  wing  configurations. 3~'  The  effect  on  aeroelastic  stability  of  such 
parameters  as  inboard  wing  folding  angle  and  hinge  stiffness  were  investigated.  Conclusions 
from  these  studies  include  l)a  trend  of  increasing  flutter  dynamic  pressure  with  increasing 
inboard  wing  folding  angle  2)higher  sensitivity  of  flutter  dynamic  pressure  with  respect  to 
outboard  hinge  stiffness  as  compared  to  the  inboard  hinge  stiffness  and  3)  morphing  wing 
actuation  energy  depends  on  center  of  gravity  position,  Mach  number  and  wing  sweep. 

Lee  and  Chen8  studied  the  effect  of  hinge  free-play  on  the  flutter  and  limit  cycle  oscillation 
of  folding  wings.  They  observed  that,  for  folding  angles  between  0  and  30  degrees,  limit  cycle 
oscillation  occurs  even  when  the  flight  altitude  is  above  that  of  the  flutter  boundary.  Further 
increases  of  the  folding  wing  angle  above  30  degrees  resulted  in  no  limit  cycle  oscillation. 
They  observed  that  when  limit  cycle  oscillation  did  occur,  it  could  occur  even  for  very  small 
values  of  the  free-play  parameter  (±0.02°). 

In  recent  work,  Tang  and  Dowell9  introduced  the  concept  of  using  component  modal  anal¬ 
ysis  to  model,  efficiently  and  accurately,  multi-component  folding  wing  configurations.  Their 
results  showed  that  an  increase  in  inboard  wing  hinge  stiffness  lead  to  improved  aeroelastic 
stability  while  an  increase  in  outboard  hinge  stiffness  decreased  the  flutter  velocity.  They 
also  noted  that  for  certain  values  of  inboard  hinge  stiffness  a  hump-type  flutter  could  occur. 
Their  computational  results  were  also  compared  to  experiment  and  in  the  experiment  limit 
cycle  oscillation  was  observed.  Good  agreement  was  found  between  computed  and  measured 
flutter  velocities  and  frequencies.  However  due  to  the  linear  nature  of  the  structural  and 
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aerodynamic  computational  models  used  in  this  work,  the  prediction  of  these  limit  cycles 
was  not  possible. 

The  work  presented  in  this  paper  is  an  extension  of  that  presented  in  Ref.  9  to  include 
nonlinear  effects  in  both  the  structural  and  aerodynamic  models.  Geometric  nonlinearity  in 
the  structure  is  modeled  using  von  Karman  strains  with  Kirchhoff  thin-plate  theory  and  the 
resulting  nonlinear  variational  statement  is  discretized  with  a  discrete  Ritz  basis  computed 
using  a  combined  finite-element/component  synthesis  analysis.  The  flow  is  modeled  using 
a  vortex  lattice  potential  model  which  accounts  for  the  nonlinear  tangent  flow  boundary 
conditions.  Post-flutter  limit  cycle  results  from  the  computational  model  are  compared  to 
those  measured  in  experiment  for  three  different  outboard  wing  folding  angle  configurations. 


II.  Folding  Wing  Configuration 


A  schematic  of  the  folding  wing  geometry  which  is  to  be  investigated  is  shown  in  Figure  1. 
The  folding  wing  system  consists  of  three  separate  components,  component  A,  i.e.  fuselage; 
component  B  ,  i.e.  inboard  wing  and  component  C  ,  i.e.  outboard  wing.  Components  A 
and  B  are  attached  through  a  hinge  which  is  modeled  as  a  set  of  torsional  springs  at  several 
points,  Pj  with  spring  stiffness,  Kaj  ■  Components  B  and  C  are  also  attached  through  a 
hinge  which  is  modeled  as  a  set  of  torsional  springs  at  several  points,  Sj  with  torsional 
spring  stiffness,  KBy  The  hinge  model  is  assumed  to  have  negligible  mass  compared  to  the 
wing  structure  model.  The  initial  folding  angles  between  components  A  and  B,  and  B  and 
C  are  9 b  and  9C,  respectively.  Note  that  9b  and  9C  are  the  static  equilibrium  angles.  These 
depend  upon  the  initial  unsprung  wing  folding  angle  and  sprung  deformation  due  to  wing 
gravity.  In  addition  to  the  main  coordinate  system,  xyz,  two  additional  local  coordinate 
systems,  CbVb^b  and  C,c9dc  for  components  B  and  C  are  used  as  shown  in  Figure  1. 

The  relationships  between  the  main  coordinate  system  and  local  coordinate  systems  are 
expressed  as  follows. 
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Figure  1.  Schematic  of  a  folding  wing  geometry  and  coordinate  systems. 

III.  Structural  Equations  of  Motion 

A.  Nonlinear  Variational  Statement 

The  statement  of  virtual  work  requires  that  the  work  done  by  externally  applied  loads  be 
equal  to  the  work  done  by  inertial,  dissipative  and  internal  forces  for  any  virtual  displacement. 
For  a  structure  with  volume  V  and  surface  area  S  this  can  be  written  as 

/  5uTFdV  +  f  Su7 &trdS  =  f  (<furpu  +  5u7  cu  +  Se1  cr )  dV  ,  (2) 

Jv  Js  Jv 

where  F  and  d>tr  represent  prescribed  body  forces  and  surface  tractions,  p  is  the  mass  density 
and  c  is  a  damping  parameter.  Also  <5u,  5e  and  cr  represent  the  vectors  of  virtual  displace¬ 
ments,  virtual  Green- Lagrange  strains  and  second  Kirchhoff-Piola  stresses  respectively. 

Here  it  will  be  assumed  that  Kirchhoff  plate  theory  applies  and  that  plate  rotations  are 
negligible  compared  to  unity  (von  Karman  assumption).  With  these  simplifications,  the 
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Green- Lagrange  strains  can  be  written  as 


1  .2 

£ XX  ^ X  “I-  2 

(3) 

1  .2 

eyy  ~  vy  +  c^y  zwyy 

(4) 

^xy  ^ y  “l-  ^ x  ^x^Uy  ‘^Z'^xy  5 

(5) 

where  to  the  degree  of  accuracy  given  by  the  von  Karman  and  Kirchhoff  assumptions,  the 
x,  y  and  z  displacements  u,v  and  w,  which  make  up  the  vector  u  in  Eq.  2,  are  given  as: 


„  dw 

(6) 

u  ~  u  —  z—— 
ox 

^  dw 

(7) 

V  PS  v  —  z—— 
dy 

w  ~  w  . 

(8) 

The  displacement  vector  u  now  can  be  considered  to  contain  only  the  mid-plane  displace¬ 
ments  u,  v  and  w. 

To  be  consistent  with  the  Kirchhoff  and  von  Karman  assumptions  it  is  assumed  that  the 
constitutive  model  is  linear,  plane  stress.  Hence,  given  a  constant  component  thickness  h, 
we  can  define  the  membrane  (N)  and  bending  (M)  resultant  forces  as: 


j-h/2 

rh/2 

N  = 

/  crdz  = 

/  Cedz  = 

J-h/2 

J-h/2 

j-h/2 

j-h/2 

M  = 

/  azdz  = 

-  /  Cezdz 

J-h/2 

J-h/2 

f/i/2 


I -h/2 


C  (em  +  zx)  dz  =  hCer 


fh/ 2  h3 

/  C  (em  +  zx)  zdz  =  —Cx  , 
J-h/2 


(9) 

(10) 


where  cr  =  Ce  is  the  linear  plane  stress  constitutive  relation  and  em  and  x  are  the  membrane 
strains  and  curvatures  respectively  (see  Eqs.  3-5). 

Integrating  the  volume  integrals  in  Eq.  2  over  the  thickness  and  substituting  Eqs.  9,  10 
gives  the  following  virtual  work  expression: 


h  (  Su1  F dS  +  f  SuT$>trdS 
Js  Js 


h  /  5uT piidS  +  h  5uTcudS  +  /  Se‘nNdS  +  /  5xTMdS  . 

Js  Js  Js  Js 

(11) 


At  this  point  a  decision  must  be  made  as  to  how  Eq.  11  is  to  be  semi-discretized  in 
space.  In  other  words,  what  sort  of  trial  functions  (Ritz  basis)  will  be  used  to  represent  the 
displacements  contained  in  the  displacement  vector  u  and  the  variation  of  this  vector  5u. 
Popular  choices  include  hat  functions  (in  finite  element  methods)  and  eigenfunctions  of  the 
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linearized  system.  Here  we  will  choose  these  functions  to  be  N-dimensional  vectors  derived 
using  component  synthesis  analysis  in  combination  with  finite  element  analysis  (to  be  more 
fully  explained  in  the  next  section). 

The  displacement  vector  D  corresponding  to  the  x,y  and  z  mid-plane  displacements  at  N 
discrete  points  in  space  and  the  variation  of  this  displacement  vector  c)D  can  now  be  written 
respectively  as 


D  =  *<(t)  , 


(12) 


and 


6 D  =  ,  (13) 

where  is  a  3N  x  P  matrix  containing  the  P  37V-dimensional  basis  vectors  and  £  is  a 
vector  of  independent  generalized  time  dependent  degrees  of  freedom.  Note,  if  instead  of 
finite  element  analysis,  continuous  functions  were  used  in  the  component  synthesis  procedure, 
N  would  be  equal  to  one. 

Each  basis  vector  in  \Ef  satisfies  the  system  geometric  boundary  conditions  (in  the  discrete 
sense)  but  not  necessarily  the  natural  boundary  conditions  and  hence  these  vectors  can  be 
considered  to  be  in  the  admissible  class.10 

The  vectors  of  membrane  strains  and  curvatures  and  virtual  membrane  strains  and  cur¬ 
vatures  are  expressed  in  terms  of  the  generalized  displacements  and  virtual  displacements  as 
(using  standard  “B  matrix”  notation11): 


£,„  =  (b;,.  +  b:‘)  c 

(14) 

Sem  =  (B(„  +  B"1)  SC 

(15) 

X  — 

(16) 

Sx  =  Bw5<  . 

(17) 

Substitution  of  Eqs.  13-17  into  Eq.  11  gives: 


5C  (  h  V1  F dS  +  /  V1  <f>trdS  -  h  /  V1  VpdS C  -h  V1  cVdS C 

Js  Js  Js  Js 


(BJ'N dS  -  I  (BZYNdS-  j  B,(M</,S')  =0  . 


(18) 


Considering  the  arbitrariness  of  the  virtual  generalized  displacements,  Eq.  18  can  be  written 
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as: 


h  I  ^T^pdSC  +  h  I  ^cVdSC  +  [  (BlJTN dS  +  [  (B”(  )T  NILS' + 

Js  Js  Js  Js 

I  B^M dS  =  h  I  VTFdS+  I  *T$trdS  .  (19) 

Js  Js  Js 

Equation  19  is  the  semi-discretized  (in  space)  equation  of  motion  for  the  system.  The 
integrals  over  dS  in  Eq.  19  are  performed  for  each  component  using  the  midpoint  rule  in  the 
local  coordinate  system  of  the  component.  In  order  to  compute  the  derivatives  needed  to 
construct  the  “B”  matrices  in  Eq.  19,  third-order  finite  difference  approximations  are  used. 

Discretization  in  time  of  Eq.  19  is  accomplished  through  the  use  of  the  Hilbert-Hughes- 
Taylor  (HHT)  implicit,  second  order,  finite  difference  scheme.11 

B.  Derivation  of  Ritz  Basis  Through  Component  Synthesis 

As  the  details  of  component  modal  synthesis  are  presented  in  great  detail  in  many  references 
(see  Refs.  9, 10, 12-18), the  general  homogeneous  disjoint  set  of  equations  for  M  components 
will  be  given  without  derivation.  This  set  of  disjoint  equations  is  generated  by  writing 
Lagrange’s  equations  for  each  of  the  components  and  when  put  into  one  matrix  system  can 
be  expressed  as 


+  C dCd  +  K-dCd  —  0  >  (20) 

where  =  [Ci  >  is  the  disjoint  generalized  displacement  vector  of  dimension  R  and 

and  are  the  R  x  R  block-diagonal  disjoint  mass,  damping  and  stiffness  matri¬ 
ces  respectively.  The  word  disjoint  is  used  here  to  describe  a  set  of  coordinates  which  are 
not  independent.  The  value  of  R  corresponds  to  the  sum  of  the  number  of  basis  vectors 
chosen  for  each  component.  After  the  component  (vector)  bases  are  generated  using  finite 
element  discretization,  the  disjoint  matrices  in  Eq.  20  are  computed  using  numerical  inte¬ 
gration  (midpoint  rule)  over  each  component  area  (for  the  in-plane  stiffness  matrices)  and 
orthonormality  (for  the  out-of-plane  stiffness  matrices).  Inertial  contributions  are  considered 
for  each  component  of  the  displacement. 

In  the  current  work,  the  potential  energy  expressions  used  to  derive  Eqs.  20  include 
contributions  from  the  linear  plate  bending  and  membrane  strain  energy  and  the  energy 
due  to  the  torsional  springs  at  the  attachment  points.  The  potential  energy  U  due  to  these 
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springs,  at  the  interfaces  A  —  B  and  B  —  C  (see  Fig.  1),  is  given  by 


v*-b  -  i  E,  JCwflfrU.  -  S&U-.)’  •  i-i . np 

=  ,  3  =  1 . NS, 

where  Np  and  Ns  are  the  total  numbers  of  the  attachment  points  of  the  two  interfaces. 

The  individual  component  generalized  displacement  vectors  Ci  can  in  turn  be  written  as 

cf  =[af,bf,qf]  ,  (21) 

where  a*.,  b,  and  b,  are  the  ith  component  generalized  coordinates  corresponding  to  the 
vectors  of  x,y  and  z  mid-plane  displacements  respectively: 

*‘  =  E<#S 

3 

3 

3 

In  Eqs.  22-24,  and  7*  correspond  to  the  jth  basis  vector  used  to  expand  the  ith 

component  x,y  and  z  mid-plane  displacement  vectors. 

In  this  paper  the  S  =  3 Np  +  3 Ns  interface  constraints  are  enforced  through  the  writing 
of  a  matrix  equation  of  dimension  S  x  R  which  is  written  as 

AO  =  0  .  (25) 

In  order  to  remove  the  redundant  coordinates  from  a  coordinate  transformation  between 
Cd  and  the  set  of  independent  coordinates  £  is  performed.  This  transformation  is  given  by 
the  equation 


(22) 

(23) 

(24) 


C  =  T<  .  (26) 

This  coordinate  transformation  consists  of  finding  a  basis  for  the  null  space  of  the  matrix 
A  and  then  expanding  C,d  in  these  vectors.  The  basis  is  computed  using  a  singular  value 
decomposition  of  the  matrix  A.  It  should  be  noted  that  an  alternative  method  of  handling 
the  constraints  would  be  through  the  introduction  of  Lagrange  multipliers  into  Lagrange’s 
equations.12 

In  the  above  expression  for  the  disjoint  equations  of  motion,  the  coordinate  transfor- 
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mation  between  the  3N  dimensional  vector  of  physical  displacements  D  and  Cd  is  given 

by 


D  =  SIC,  ,  (27) 

where  Q,  is  a  matrix  composed  of  all  the  component  basis  vectors.  Using  the  coordinate 
transformation  given  in  Eq.  26,  Eq.  27  can  be  written  as 

D  =  f2TC  =  T'C  •  (28) 

The  matrix  of  global  basis  vectors  \l/  is  then  used  in  the  nonlinear  variational  statement, 
Eq.  19. 

A  further  reduction  in  dimension  can  be  performed  if  Eq.  26  is  placed  into  Eq.  20  and 
the  system  is  then  premultiplied  by  the  transpose  of  T.19  The  resulting  system  can  then  be 
cast  as  an  eigevalue  problem  and  the  matrix  of  eigenmodes  T  then  used  in  the  coordinate 
transformation: 


D  =  QTTCr  =  VrCr  ■  (29) 

While  this  procedure  can  lead  to  a  drastically  reduced  system  dimension,  unfortunately 
the  resulting  global  basis  vectors  satisfy  (in  the  discrete  sense)  the  linear  natural  in-plane 
boundary  conditions  which  degrades  the  solution  convergence  of  the  nonlinear  problem  with 
respect  to  the  required  number  of  basis  functions.20 

C.  Component  Basis  Selection 

The  components  of  the  folding  wing  structure  studied  here  have  the  following  material/geometric 
properties.  Component  A  is  constructed  from  plcxiglas  with  a  modulus  3.05  x  109  Pa.,  den¬ 
sity  1145  kg/m3,  Poisson’s  ratio  of  0.45  and  thickness  0.00238  m.  Components  B  and  C  are 
constructed  from  aluminum  with  modulus  7.2  x  1010  Pa.,  density  2850  kg/m3,  Poisson’s  ratio 
0.3  and  thickness  0.000254  m.  The  stiffness  values  Kaj,  Kbj  used  for  the  torsional  springs  at 
the  component  interfaces  is  0.18  kg-rn/s2.  For  the  results  presented  in  this  work,  the  value 
of  9b  is  30  degrees  and  three  values  of  9C  are  investigated,  0,  30  and  60  degrees. 

The  interface  constraints  (as  expressed  in  Eq.  25)  are  that  the  component  mid-plane 
deflections  ( u ,  v,  w)  in  the  global  xyz  coordinate  system  are  to  be  the  same  at  the  hinge 
attachment  points. 

Each  of  the  component  bases  used  in  this  work  are  computed  using  a  finite  element  modal 
analysis.  The  commercial  code  ANSYS21  is  used  for  this  purpose.  The  following  types  of 
component  bases  are  used  to  expand  the  out-of-plane  deflections  w: 
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Component  A:  the  clastic  out-of-plane  clamped-free-free-free  modes  from  a  modal  analysis 
of  Component  A. 

Component  B:  the  out-of-plane  elastic  and  rigid  body  modes  from  a  free- free- free- free 
modal  analysis  of  Component  B. 

Component  C:  the  out-of-plane  elastic  and  rigid  body  modes  from  a  free-free-free-free 
modal  analysis  of  Component  C. 

The  component  in-plane  deflections  are  expanded  using  the  following  types: 

Component  A:  both  u  and  v  are  expanded  using  the  elastic  out-of-plane  simply-supported- 
free-free-free  modes  computed  with  a  modal  analysis  of  component  A. 

Component  B:  both  u  and  v  are  expanded  using  the  3  rigid  body  modes  from  a  linear 
plane  stress  modal  analysis  of  component  B  and  the  remaining  basis  vectors  are  the  clastic 
out-of-plane  free-free-free-free  modes  for  component  B 

Component  C:  both  u  and  v  are  expanded  using  the  3  rigid  body  modes  from  a  linear 
plane  stress  modal  analysis  of  component  C  and  the  remaining  basis  vectors  are  the  clastic 
out-of-plane  free-free-free-free  modes 

The  number  of  out-of-plane  modes  used  for  components  A,  B  and  C  is  25,  75  and  75 
respectively  while  the  number  of  in-plane  modes  is  5,  15  and  15.  These  choices  ,  after 
the  redundant  coordinates  are  eliminated,  give  a  total  number  of  structural  degrees  of  free¬ 
dom  of  221.  A  convergence  study  showed  that  this  number  of  modes  gave  a  good  solution 
convergence. 

It  should  be  noted  that  the  choice  of  the  types  of  bases  used  in  the  expansion  of  the 
in-plane  components  is  guided  by  the  work  presented  in  Ref.  20. 

The  results  for  the  first  five  modal  frequencies  for  the  three  outboard  wing  folding  config¬ 
urations  are  shown  in  Table  1  (denoted  as  Theory  in  Table  1)  along  with  the  corresponding 
results  from  an  ANSYS  modal  analysis  conducted  using  a  finite  element  discretization  of  the 
complete  structure  and  experiment. 

IV.  Aerodynamic  Equations 

The  flow  held  about  the  folding  wing  model  is  assumed  to  be  potential  how.  Therefore 
the  equations  of  motion  for  the  fluid  can  be  reduced  to  Laplace’s  equation  for  the  velocity 
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ec  =  o 

9C  =  30 

9C  =  60 

Mode 

Present 

Theory 

ANSYS 

Test 

Present 

Theory 

ANSYS 

Test 

Present 

Theory 

ANSYS 

Test 

1 

4.67 

4.66 

4.75 

4.53 

4.54 

4.25 

4.67 

4.67 

4.0 

2 

16.67 

16.82 

16.0 

17.01 

17.18 

17.25 

16.67 

16.77 

15.25 

3 

22.23 

22.88 

32.5 

36.40 

37.04 

31.50 

22.22 

22.50 

23.50 

4 

41.72 

42.15 

43.75 

43.39 

43.63 

46.25 

41.71 

41.80 

47.5 

5 

66.17 

67.68 

75.25 

67.99 

69.45 

73.25 

66.15 

67.19 

75.25 

Table  1.  Comparison  of  first  five  folding  wing  modal  frequencies. 


potential  $ 


V2$  =  0  .  (30) 

The  boundary  conditions  which  must  be  satisfied  are  those  of  zero  normal  flow  on  the  wing,i.e. 

(V<&  +  q)  •  n  =  0  ,  (31) 

and  also  that  the  disturbance  created  by  the  potential  must  decay  at  distances  far  from  the 
wing.  The  latter  equation  can  be  expressed  as 

lim  V4>  =  0  .  (32) 

r— xx) 

In  Eq.  31  q  is  the  relative  velocity  between  the  undisturbed  fluid  in  the  fluid  domain  and  the 
wing.  Using  Green’s  identity  it  can  be  shown22  that  a  solution  to  Eqs.  30-32  can  be  found  by 
distributing  elementary  solutions  to  Laplace’s  equation  on  the  problem  boundaries.  For  the 
model  presented  here  this  is  accomplished  by  distributing  vortex  rings  on  the  wing  surface 
and  in  the  wake. 

In  this  work  the  small  disturbance  assumption  is  not  made  and  hence  the  aerodynamic 
grid  used  to  define  bound  collocation  and  vortex  locations  is  moved  in  three  dimensions  using 
the  three  components  of  the  displacment  vector.  It  was  determined  through  preliminary 
calculations  that  imposing  the  physically  correct,  but  numerically  expensive  force-free  wake 
condition  had  very  little  effect  on  the  solution  (see  Fig.  2)  hence  the  results  presented  here 
use  a  planar  wake  assumption.  However  the  trailing  edge  wake  shedding  location  is  taken 
to  be  the  deformed  position.  See  Katz  and  Plotkin23  for  further  details  on  the  use  of  vortex 
rings  to  discretize  the  flow  model. 

The  pressure  on  the  wing  p,  which  is  the  external  load  acting  on  the  structure  and  used  to 
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compute  &tr 
equation: 


in  Eq.  19,  is  calculated  using  the  unsteady  (and  nonlinear)  form  of  Bernoulli’s 


Poo  -p 
Poo 


(V4>)2 

2 


—  q  •  V<h  + 


~dt 


(33) 


where  and  are  the  free  stream  pressure  and  density  respectively.  The  time  derivative 
in  Eq.  33  is  computed  using  a  first-order  backward  difference  approximation. 


V.  Coupling  of  Structural  and  Fluid  Models 

Coupling  of  the  aerodynamics  to  the  structural  dynamics  occurs  through  the  specification 
of  the  tractions  in  Eq.  19  and  by  the  resulting  structural  deflection  and  velocity  ,  which 
is  used  in  the  solution  of  Eqs.  30,  31. 

Strong  coupling  of  the  fluid  and  structural  models  is  achieved  via  subiteration  within 
each  timestep  of  the  simulation.24  As  the  grid  used  in  the  computation  of  the  component 
bases  is  not  coincident  with  the  grid  used  in  defining  the  aerodynamic  collocation  points, 
the  tractions,  structural  displacements  and  structural  velocities  are  interpolated  using  local 
bilinear  interpolation.25 


VI.  Theoretical  and  Experimental  Correlations 

The  results  presented  in  this  section  are  computed  using  the  following  numerical  pa¬ 
rameters.  The  number  of  chordwise  aerodynamic  bound  vortex  rings  is  30  and  the  number 
spanwise  bound  vortex  rings  is  33.  A  total  of  114  panels  are  used  in  the  chordwise  direction 
for  the  wake.  With  a  timestep  of  0.0005  seconds  and  a  minimum  flow  velocity  investigated 
of  15.00  m/s,  this  corresponds  to  a  minimum  of  3  chord  lengths  of  wake.  Structural  damping 
is  not  used  in  the  simulation,  however  numerical  damping  is  present  in  the  HHT  scheme. 
The  damping  parameter  in  the  HHT  scheme  is  chosen  to  be  0.05  which  has  been  found  to 
give  good  results  for  nonlinear  problems.26,27  All  response  results  shown  in  this  section  are 
for  the  structural  deflection/velocity  in  the  global  z  coordinate  direction. 

Figures  2-4  compare  the  inboard  wing  (Component  B),  trailing  edge  tip  root-mean  square 
velocity  of  the  limit  cycle  oscillations  (LCO)  from  the  current  computational  model  to  those 
measured  in  the  wind-tunnel  experiments.  Results  for  three  different  folding  wing  configura¬ 
tions  are  presented.  For  each  of  the  configurations,  the  correlation  between  computed  and 
measured  resonse  at  this  point  on  the  folding  wing  is  good.  It  does  appear,  however,  that 
the  computational  model  slightly  underpredicts  the  point  of  bifurcation  (“flutter  velocity”). 
The  computational  model  predicts  flutter  velocities  of  15.5  m/s, 15. 35  m/s  and  15.8  m/s  for 
the  configurations  corresponding  to  9C  =  0 ,9C  =  30  and  9C  =  60  respectively.  While  difficult 
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Figure  2.  Component  B  trailing  edge, tip  root-mean  square  LCO  velocity  versus  flow  velocity; 

eb  =  30, 9C  =  0. 


to  determine  exactly,  the  flutter  velocities  in  the  experiment  for  these  configurations  are 
approximately  18  m/s,  17. 5  m/s  and  17  m/s. 

Also  shown  in  Fig.  2  is  a  computation,  for  a  flow  velocity  of  16  m/s,  which  includes  an 
aerodynamic  model  which  uses  both  the  exact  nonlinear  tangent  flow  boundary  condition 
and  the  imposition  of  the  force- wake  free  condition  (denoted  as  THEORY:  NL  STRUC, 
NL  BC,  FF  WAKE  in  the  figure)  and  a  computation  at  16  m/s  which  uses  the  linearized 
tangent  flow  boundary  condition  and  a  planar  wake(denoted  as  THEORY:  NL  STRUC,  LIN 
BC,  PLANAR  WAKE  in  the  figure).  While  a  definite  conclusion  can  only  be  drawn  for 
the  flow  velocity  of  16  m/s  and  6b  =  30, 9C  =  0  folwing  wing  configuration,  it  appears  that 
for  the  stated  flow  velocity  and  folding  wing  configuration, including  the  nonlinear  boundary 
condition  in  the  model  appears  to  have  a  slight  effect  while  the  inclusion  of  the  force-free  wake 
model  has  minimal  effect.  Results  presented  in  the  remainder  of  this  section  are  computed 
using  a  model  which  includes  a  nonlinear  structural  model,  the  exact  nonlinear  tangent  flow 
boundary  condition  and  a  planar  wake  assumption. 

A  comparison  of  the  theoretical  and  experimental  LCO  dominant  response  frequencies  as 
a  function  of  the  flow  velocity,  as  computed  from  a  discrete  Fourier  transform  of  the  response 
time  history,  is  shown  in  Figs.  5-7.  For  each  of  the  configurations,  near  the  flutter  velocity 
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Figure  3.  Component  B  trailing  edge, tip  root-mean  square  LCO  velocity  versus  flow  velocity; 

eb  =  30, ec  =  30. 


Figure  4.  Component  B  trailing  edge, tip  root-mean  square  LCO  velocity  versus  flow  velocity; 

eb  =  30, 9C  =  60. 
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FLOW  VELOCITY  (M/S) 

Figure  5.  LCO  dominant  response  frequency  for  component  B  trailing  edge  tip  velocity  versus 
flow  velocity;  8 j,  =  30, 8C  =  0. 


FLOW  VELOCITY  (M/S) 

Figure  6.  LCO  dominant  response  frequency  for  component  B  trailing  edge  tip  velocity  versus 
flow  velocity;  6b  =  30, 8C  =  30. 
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Figure  7.  LCO  dominant  response  frequency  for  component  B  trailing  edge  tip  velocity  versus 
flow  velocity;  6 &  =  30, 9C  =  60. 


the  comparison  between  the  computation  and  experiment  for  these  response  frequencies  is 
good.  Also,  in  both  the  computation  and  experiment,  the  response  frequencies  appear  to 
increase  with  increasing  flow  velocity.  However,  unlike  the  experiment,  as  the  flow  velocity 
increases,  the  computational  model  results  show  increasing  dynamic  complexity  in  the  form 
of  multiple  frequencies  having  significant  contributions  to  the  response.  This  behavior  is 
shown  in  Fig.  8  which  compares  the  discrete  Fourier  transform  for  the  computational  and 
experimental  model  at  a  flow  velocity  of  19.723  m/s  for  the  the  9C  =  0  configuration.  The 
computational  results  show  the  presence  of  multiple  contributing  response  frequencies, many 
of  which  are  not  integer  multiples  of  each  other.  This  increase  in  dynamic  complexity  causes 
the  dominant  frequency  to  switch  to  a  higher  frequency  for  higher  flow  velocities  as  is  seen 
in  Figs.  5  and  6.  The  complex  nature  of  the  computational  results  is  further  shown  in  Figs.  9 
and  10  which  give  the  velocity  response  time  history  and  phase  plots  for  the  component 
B  tip, trailing  edge  location  at  a  flow  velocity  of  22.70  m/s.  The  folding  wing  configuration 
for  these  two  figures  is  9C  =  0.  In  contrast  to  Fig.  10,  Fig.  11  is  a  phase  plot  for  the  same 
configuration  at  a  flow  velocity  of  16  m/s  which  is  near  the  flutter  velocity.  The  differences 
in  the  nature  of  the  phase  portraits  are  apparent. 

In  order  to  attempt  to  quantify  this  dynamic  complexity  a  measure  has  been  used  which 
is  computed  as  the  ratio  of  the  magnitude  of  the  maximum  Fourier  coefficient  Amax  to  the 
magnitude  of  all  of  the  Fourier  coefficients  where  this  latter  value  is  computed  as  y/JA  Af. 
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FREQUENCY  (Hz) 


Figure  8.  Discrete  Fourier  transform  of  component  B  trailing  edge  tip  LCO  velocity  at  a  flow 
velocity  of  19.723  m/s;  0C  =  0. 


Figure  9.  Component  B  trailing  edge, tip  LCO  velocity  versus  time  at  a  flow  velocity  of  22.70 

m/s;  6C  =  0. 
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In  order  to  account  for  discrete  sampling  the  maximum  and  the  two  surrounding  points  are 
used  in  the  computation  of  Amax.  A  small  value  of  this  measure  denotes  higher  dynamic 
complexity  in  the  form  of  more  freauencv  components  contributing:  significantly  to  the  resop- 
nse.  T  L4. 

As  cj  )w 


Figure  12.  Component  B  tip, trailing  edge  velocity  dynamic  complexity  measure  versus  flow 
velocity;  6/,  =  30, 8C  =  0. 


velocity  while  the  experimental  models  appear  to  have  mostly  a  single  frequency  response. 

It  is  interesting  that  the  outboard  wing  computational  results  shown  in  Fig.  15  do  not 
show  this  trend.  This  figure  demonstrates  that  for  the  outboard  wing  a  single  frequency 
response  occurs  over  the  flow  velocity  range  investigated.  Also  from  Fig.  16  it  can  be  seen 
that  the  computed  tip  displacement  of  the  outboard  wing  is  quite  large,  especially  for  the 
two  configurations  with  the  smaller  outboard  wing  folding  angle  (0C  =  0  and  9C  =  30).  The 
thetac  =  60  configuration  shows  much  smaller  outboard  wing  tip  limit  cycle  amplitudes  and 
a  larger  degree  of  “nonlinear  stiffening”  as  the  flow  velocity  is  increased. 

A.  Conclusions 

The  nonlinear  aeroelastic  response  characteristics  of  a  three  component  (fuselage,  inboard 
wing  and  outboard  wing)  folding  wing  configuration  are  investigated  both  experimentally 
and  theoretically.  The  theoretical  analysis  includes  a  nonlinear  structural  dynamics  model 
using  a  discrete  Ritz  basis  derived  from  finite  clmcnt  and  component  synthesis  analysis  and 
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Figure  13.  Component  B  tip, trailing  edge  velocity  dynamic  complexity  measure  versus  flow 
velocity;  6b  =  30, 8C  =  30. 


'  10  15  20  25 
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Figure  14.  Component  B  tip, trailing  edge  velocity  dynamic  complexity  measure  versus  flow 
velocity;  6b  =  30, 8C  =  60. 
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Figure  15.  Theoretical  component  C  tip,  trailing  edge  velocity  dynamic  complexity  measure 
versus  flow  velocity. 


Figure  16.  Theoretical  component  C  tip,  trailing  edge  root  mean  square  LCO  displacement 
versus  flow  velocity. 
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a  nonlinear  vortex  lattice  potential  flow  model  for  the  fluid. 

Computational  results  for  an  inboard  wing  folding  angle  of  30  degrees  and  three  different 
outboard  wing  folding  angles  (0,30  and  60  degrees)  are  compared  to  experiment.  For  each 
configuration,  post-flutter  limit  cycle  oscillations  are  found  in  both  the  experiment  and 
computations.  Overall  the  correlation  between  the  computation  and  experiment  in  this 
post-flutter  region  is  good.  Predicted  and  experimental  response  amplitude  and  frequencies 
for  the  inboard  wing  trailing  edge  tip  position  compare  favorably.  The  theoretical  model 
predicts  similar  outboard  limit  cycle  tip  displacements  for  outboard  wing  folding  angles  of  0 
and  30  degrees  while  the  configuration  where  the  outboard  folding  angle  is  60  degrees  show 
signficantly  smaller  limit  cycle  amplitudes  and  a  higher  degree  of  nonlinear  stiffening  with 
increasing  flow  velocity. 

A  measure  of  dynamic  complexity  consisting  of  the  ratio  of  the  maximum  Fourier  co¬ 
efficient  amplitude  to  the  square  root  of  the  sum  of  the  squares  of  all  the  amplitudes  is 
introduced.  From  this  measure,  and  by  inspection  of  the  discrete  Fourier  transforms,  time 
histories  and  phase  portaits,  it  is  found  that  the  computational  model  inboard  wing  response 
contains  larger  contributions  from  multiple  frequencies  as  compared  to  the  experiment  lead¬ 
ing  to  a  more  complex  dynamic  response.  The  omission  of  structural  dynamic  damping  or 
an  overly  simplified  theoretical  model  of  the  joint  dynamics  is  a  possible  explain  for  this 
discrepancy.  This  complex  dynamic  response  is  not  predicted  for  the  outboard  wing,  as  the 
limit  cycle  behavior  consists  mostly  of  a  single  frequency,  large  amplitude,  response. 

Finally  a  significant  result  from  this  study  is  the  determination  that  both  the  experiment 
and  computation  show  that  differences  in  the  limit  cycle  behavior  do  exist  between  the 
various  folding  wing  angles.  Also,  based  upon  the  calculations  using  a  linear  aerodynamic 
model,  it  appears  that  the  nonlinear  aerodynamic  effects  are  smaller  than  the  nonlinear 
structural  effects  for  the  cases  studied. 
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